home *** CD-ROM | disk | FTP | other *** search
/ Info-Mac 11 / Info-Mac_XI_Disc_1.cdr_ / Info-Mac XI Disc 1.cdr / Programs / Science & Math / MacEspresso 1.0 / espresso / mincov.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-02-26  |  9.6 KB  |  370 lines  |  [TEXT/R*ch]
  1. #include "mincov_int.h"
  2.  
  3. /*
  4.  *  mincov.c
  5.  */
  6.  
  7. #define USE_GIMPEL
  8. #define USE_INDEP_SET
  9.  
  10. extern int select_column();
  11. extern void select_essential();
  12. extern int verify_cover();
  13.  
  14. #define fail(why) {\
  15.     (void) fprintf(stderr, "Fatal error: file %s, line %d\n%s\n",\
  16.     __FILE__, __LINE__, why);\
  17.     (void) fflush(stdout);\
  18.     abort();\
  19. }
  20.  
  21. sm_row *
  22. sm_minimum_cover(A, weight, heuristic, debug_level)
  23. sm_matrix *A;
  24. int *weight;
  25. int heuristic;        /* set to 1 for a heuristic covering */
  26. int debug_level;    /* how deep in the recursion to provide info */
  27. {
  28.     stats_t stats;
  29.     solution_t *best, *select;
  30.     sm_row *prow, *sol;
  31.     sm_col *pcol;
  32.     sm_matrix *dup_A;
  33.     int nelem, bound;
  34.     double sparsity;
  35.  
  36.     /* Avoid sillyness */
  37.     if (A->nrows <= 0) {
  38.     return sm_row_alloc();        /* easy to cover */
  39.     }
  40.  
  41.     /* Initialize debugging structure */
  42.     stats.start_time = util_cpu_time();
  43.     stats.debug = debug_level > 0;
  44.     stats.max_print_depth = debug_level;
  45.     stats.max_depth = -1;
  46.     stats.nodes = 0;
  47.     stats.component = stats.comp_count = 0;
  48.     stats.gimpel = stats.gimpel_count = 0;
  49.     stats.no_branching = heuristic != 0;
  50.     stats.lower_bound = -1;
  51.  
  52.     /* Check the matrix sparsity */
  53.     nelem = 0;
  54.     sm_foreach_row(A, prow) {
  55.     nelem += prow->length;
  56.     }
  57.     sparsity = (double) nelem / (double) (A->nrows * A->ncols);
  58.  
  59.     /* Determine an upper bound on the solution */
  60.     bound = 1;
  61.     sm_foreach_col(A, pcol) {
  62.     bound += WEIGHT(weight, pcol->col_num);
  63.     }
  64.  
  65.     /* Perform the covering */
  66.     select = solution_alloc();
  67.     dup_A = sm_dup(A);
  68.     best = sm_mincov(dup_A, select, weight, 0, bound, 0, &stats);
  69.     sm_free(dup_A);
  70.     solution_free(select);
  71.  
  72.     if (stats.debug) {
  73.     if (stats.no_branching) {
  74.         (void) printf("**** heuristic covering ...\n");
  75.         (void) printf("lower bound = %d\n", stats.lower_bound);
  76.     }
  77.     (void) printf("matrix     = %d by %d with %d elements (%4.3f%%)\n",
  78.         A->nrows, A->ncols, nelem, sparsity * 100.0);
  79.     (void) printf("cover size = %d elements\n", best->row->length);
  80.     (void) printf("cover cost = %d\n", best->cost);
  81.     (void) printf("time       = %s\n", 
  82.             util_print_time(util_cpu_time() - stats.start_time));
  83.     (void) printf("components = %d\n", stats.comp_count);
  84.     (void) printf("gimpel     = %d\n", stats.gimpel_count);
  85.     (void) printf("nodes      = %d\n", stats.nodes);
  86.     (void) printf("max_depth  = %d\n", stats.max_depth);
  87.     }
  88.  
  89.     sol = sm_row_dup(best->row);
  90.     if (! verify_cover(A, sol)) {
  91.     fail("mincov: internal error -- cover verification failed\n");
  92.     }
  93.     solution_free(best);
  94.     return sol;
  95. }
  96.  
  97. /*
  98.  *  Find the best cover for 'A' (given that 'select' already selected);
  99.  *
  100.  *    - abort search if a solution cannot be found which beats 'bound'
  101.  *
  102.  *    - if any solution meets 'lower_bound', then it is the optimum solution
  103.  *      and can be returned without further work.
  104.  */
  105.  
  106. solution_t * 
  107. sm_mincov(A, select, weight, lb, bound, depth, stats)
  108. sm_matrix *A;
  109. solution_t *select;
  110. int *weight;
  111. int lb;
  112. int bound;
  113. int depth;
  114. stats_t *stats;
  115. {
  116.     sm_matrix *A1, *A2, *L, *R;
  117.     sm_element *p;
  118.     solution_t *select1, *select2, *best, *best1, *best2, *indep;
  119.     int pick, lb_new, debug;
  120.  
  121.     /* Start out with some debugging information */
  122.     stats->nodes++;
  123.     if (depth > stats->max_depth) stats->max_depth = depth;
  124.     debug = stats->debug && (depth <= stats->max_print_depth);
  125.  
  126.     /* Apply row dominance, column dominance, and select essentials */
  127.     select_essential(A, select, weight, bound);
  128.     if (select->cost >= bound) {
  129.     return NIL(solution_t);
  130.     }
  131.  
  132.     /* See if gimpel's reduction technique applies ... */
  133. #ifdef USE_GIMPEL
  134.     if ( weight == NIL(int)) {    /* hack until we fix it */
  135.     if (gimpel_reduce(A, select, weight, lb, bound, depth, stats, &best)) {
  136.         return best;
  137.     }
  138.     }
  139. #endif
  140.  
  141. #ifdef USE_INDEP_SET
  142.     /* Determine bound from here to final solution using independent-set */
  143.     indep = sm_maximal_independent_set(A, weight);
  144.  
  145.     /* make sure the lower bound is monotonically increasing */
  146.     lb_new = MAX(select->cost + indep->cost, lb);
  147.     pick = select_column(A, weight, indep);
  148.     solution_free(indep);
  149. #else
  150.     lb_new = select->cost + (A->nrows > 0);
  151.     pick = select_column(A, weight, NIL(solution_t));
  152. #endif
  153.  
  154.     if (depth == 0) {
  155.     stats->lower_bound = lb_new + stats->gimpel;
  156.     }
  157.  
  158.     if (debug) {
  159.         (void) printf("ABSMIN[%2d]%s", depth, stats->component ? "*" : " ");
  160.         (void) printf(" %3dx%3d sel=%3d bnd=%3d lb=%3d %12s ",
  161.             A->nrows, A->ncols, select->cost + stats->gimpel, 
  162.         bound + stats->gimpel, lb_new + stats->gimpel, 
  163.         util_print_time(util_cpu_time()-stats->start_time));
  164.     }
  165.  
  166.     /* Check for bounding based on no better solution possible */
  167.     if (lb_new >= bound) {
  168.     if (debug) (void) printf("bounded\n");
  169.     best = NIL(solution_t);
  170.  
  171.  
  172.     /* Check for new best solution */
  173.     } else if (A->nrows == 0) {
  174.     best = solution_dup(select);
  175.     if (debug) (void) printf("BEST\n");
  176.     if (stats->debug && stats->component == 0) {
  177.             (void) printf("new 'best' solution %d at level %d (time is %s)\n", 
  178.         best->cost + stats->gimpel, depth, 
  179.         util_print_time(util_cpu_time() - stats->start_time));
  180.         }
  181.  
  182.  
  183.     /* Check for a partition of the problem */
  184.     } else if (sm_block_partition(A, &L, &R)) {
  185.     /* Make L the smaller problem */
  186.     if (L->ncols > R->ncols) {
  187.         A1 = L;
  188.         L = R;
  189.         R = A1;
  190.     }
  191.     if (debug) (void) printf("comp %d %d\n", L->nrows, R->nrows);
  192.     stats->comp_count++;
  193.  
  194.     /* Solve problem for L */
  195.     select1 = solution_alloc();
  196.     stats->component++;
  197.     best1 = sm_mincov(L, select1, weight, 0, 
  198.                     bound-select->cost, depth+1, stats);
  199.     stats->component--;
  200.     solution_free(select1);
  201.     sm_free(L);
  202.  
  203.     /* Add best solution to the selected set */
  204.     if (best1 == NIL(solution_t)) {
  205.         best = NIL(solution_t);
  206.     } else {
  207.         for(p = best1->row->first_col; p != 0; p = p->next_col) {
  208.         solution_add(select, weight, p->col_num);
  209.         }
  210.         solution_free(best1);
  211.  
  212.         /* recur for the remaining block */
  213.         best = sm_mincov(R, select, weight, lb_new, bound, depth+1, stats);
  214.     }
  215.     sm_free(R);
  216.  
  217.     /* We've tried as hard as possible, but now we must split and recur */
  218.     } else {
  219.     if (debug) (void) printf("pick=%d\n", pick);
  220.  
  221.         /* Assume we choose this column to be in the covering set */
  222.     A1 = sm_dup(A);
  223.     select1 = solution_dup(select);
  224.     solution_accept(select1, A1, weight, pick);
  225.         best1 = sm_mincov(A1, select1, weight, lb_new, bound, depth+1, stats);
  226.     solution_free(select1);
  227.     sm_free(A1);
  228.  
  229.     /* Update the upper bound if we found a better solution */
  230.     if (best1 != NIL(solution_t) && bound > best1->cost) {
  231.         bound = best1->cost;
  232.     }
  233.  
  234.     /* See if this is a heuristic covering (no branching) */
  235.     if (stats->no_branching) {
  236.         return best1;
  237.     }
  238.  
  239.     /* Check for reaching lower bound -- if so, don't actually branch */
  240.     if (best1 != NIL(solution_t) && best1->cost == lb_new) {
  241.         return best1;
  242.     }
  243.  
  244.         /* Now assume we cannot have that column */
  245.     A2 = sm_dup(A);
  246.     select2 = solution_dup(select);
  247.     solution_reject(select2, A2, weight, pick);
  248.         best2 = sm_mincov(A2, select2, weight, lb_new, bound, depth+1, stats);
  249.     solution_free(select2);
  250.     sm_free(A2);
  251.  
  252.     best = solution_choose_best(best1, best2);
  253.     }
  254.  
  255.     return best;
  256. }
  257.  
  258. static int 
  259. select_column(A, weight, indep)
  260. sm_matrix *A;
  261. int *weight;
  262. solution_t *indep;
  263. {
  264.     register sm_col *pcol;
  265.     register sm_row *prow, *indep_cols;
  266.     register sm_element *p, *p1;
  267.     double w, best;
  268.     int best_col;
  269.  
  270.     indep_cols = sm_row_alloc();
  271.     if (indep != NIL(solution_t)) {
  272.     /* Find which columns are in the independent sets */
  273.     for(p = indep->row->first_col; p != 0; p = p->next_col) {
  274.         prow = sm_get_row(A, p->col_num);
  275.         for(p1 = prow->first_col; p1 != 0; p1 = p1->next_col) {
  276.         (void) sm_row_insert(indep_cols, p1->col_num);
  277.         }
  278.     }
  279.     } else {
  280.     /* select out of all columns */
  281.     sm_foreach_col(A, pcol) {
  282.         (void) sm_row_insert(indep_cols, pcol->col_num);
  283.     }
  284.     }
  285.  
  286.     /* Find the best column */
  287.     best_col = -1;
  288.     best = -1;
  289.  
  290.     /* Consider only columns which are in some independent row */
  291.     sm_foreach_row_element(indep_cols, p1) {
  292.     pcol = sm_get_col(A, p1->col_num);
  293.  
  294.     /* Compute the total 'value' of all things covered by the column */
  295.     w = 0.0;
  296.     for(p = pcol->first_row; p != 0; p = p->next_row) {
  297.         prow = sm_get_row(A, p->row_num);
  298.         w += 1.0 / ((double) prow->length - 1.0);
  299.     }
  300.  
  301.     /* divide this by the relative cost of choosing this column */
  302.     w = w / (double) WEIGHT(weight, pcol->col_num);
  303.  
  304.     /* maximize this ratio */
  305.     if (w  > best) {
  306.         best_col = pcol->col_num;
  307.         best = w;
  308.     }
  309.     }
  310.  
  311.     sm_row_free(indep_cols);
  312.     return best_col;
  313. }
  314.  
  315. static void 
  316. select_essential(A, select, weight, bound)
  317. sm_matrix *A;
  318. solution_t *select;
  319. int *weight;
  320. int bound;            /* must beat this solution */
  321. {
  322.     register sm_element *p;
  323.     register sm_row *prow, *essen;
  324.     int delcols, delrows, essen_count;
  325.  
  326.     do {
  327.     /*  Check for dominated columns  */
  328.     delcols = sm_col_dominance(A, weight);
  329.  
  330.     /*  Find the rows with only 1 element (the essentials) */
  331.     essen = sm_row_alloc();
  332.     sm_foreach_row(A, prow) {
  333.         if (prow->length == 1) {
  334.         (void) sm_row_insert(essen, prow->first_col->col_num);
  335.         }
  336.     }
  337.  
  338.     /* Select all of the elements */
  339.     sm_foreach_row_element(essen, p) {
  340.         solution_accept(select, A, weight, p->col_num);
  341.         /* Make sure solution still looks good */
  342.         if (select->cost >= bound) {
  343.         sm_row_free(essen);
  344.         return;
  345.         }
  346.     }
  347.     essen_count = essen->length;
  348.     sm_row_free(essen);
  349.  
  350.     /*  Check for dominated rows  */
  351.     delrows = sm_row_dominance(A);
  352.  
  353.     } while (delcols > 0 || delrows > 0 || essen_count > 0);
  354. }
  355.  
  356. static int 
  357. verify_cover(A, cover)
  358. sm_matrix *A;
  359. sm_row *cover;
  360. {
  361.     sm_row *prow;
  362.  
  363.     sm_foreach_row(A, prow) {
  364.     if (! sm_row_intersects(prow, cover)) {
  365.         return 0;
  366.     }
  367.     }
  368.     return 1;
  369. }
  370.